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Gravitational waves from binary neutron star (BNS) and black hole/neutron star (BHNS) inspi¬ 
rals are primary sources for detection by the Advanced Laser Interferometer Gravitational-Wave 
Observatory. The tidal forces acting on the nentron stars induce changes in the phase evolution of 
the gravitational waveform, and these changes can be used to constrain the nuclear equation of state. 

Current methods of generating BNS and BHNS waveforms rely on either computationally challeng¬ 
ing full 3D hydrodynamical simulations or approximate analytic solutions. We introduce a new 
method for computing inspiral waveforms for BNS/BHNS systems by adding the post-Newtonian 
(PN) tidal effects to full numerical simulations of binary black holes (BBHs), effectively replacing 
the nontidal terms in the PN expansion with BBH results. Comparing a waveform generated with 
this method against a full hydrodynamical simulation of a BNS inspiral yields a phase difference of 
< 1 radian over ~ 15 orbits. The numerical phase accuracy required of BNS simulations to measure 
the accuracy of the method we present here is estimated as a function of the tidal deformability 
parameter A. 


I. INTRODUCTION 

In September 2015, the Advanced Laser Interferom¬ 
eter Gravitational-Wave Observatory(aLIGO) directly 
detected, for the first time ever, gravitational waves 
(GWs) [1] and the network of observatories will be joined 
shortly by advanced Virgo [2] and KAGRA [3] . The most 
likely GW sources for these detectors are mergers of bi¬ 
naries consisting of neutron stars (NSs) or black holes 
(BHs) [4]. If both objects in the binary are NSs (BNS), 
or if one is a NS and the other is a BH (a BHNS binary), 
then the tidal deformability of the NS will alter the GW 
signal in a way that is dependent upon the NS equation 
of state (EOS), allowing these observatories to constrain 
the EOS [5-13]. It is therefore of key importance to un¬ 
derstand and model the influence of tidal effects on BNS 
and BHNS waveforms. We show here that a binary black 
hole (BBH) waveform can be augmented with PN tidal 
effects to accurately model a BNS system during the in¬ 
spiral portion of the binary evolution. In principle, this 
method should also be applicable to BHNS systems. 

BNS waveforms are typically computed using post- 
Newtonian (PN) methods, which are perturbative expan¬ 
sions in the invariant velocity v = {Mdcjj/dty^^, where 
M is the total mass of the system and (j) is the orbital 
phase (here we assume c = G = 1). For binaries con¬ 
sisting of nonspinning point particles, the expansion is 
known through 3.5PN order [14]. The static NS tidal 


effects first enter at 5PN order and depend upon the 
tidal deformability [15]. The parameter A^ measures 
how much each NS i deforms in the presence of a tidal 
field, and depends on the NS mass and EOS implicitly 
through its dimensionless Love number k 2 ,i and radius 
Ri'. Xi = (2/3)fc2,ii?f [5]. As V increases throughout the 
inspiral, the missing 4PN, 4.5PN, and 5PN point-particle 
terms can result in the late portion of the PN wave¬ 
form becoming inaccurate before the static tidal terms 
are large enough to contribute. For estimating the NS 
tidal deformability by using PN waveforms, the error in¬ 
troduced by neglecting the higher order PN terms can 
be as large as the statistical errors due to noise in the 
measured signals [11, 12, 16, 17]. 

Effective-one-body (EOB) models that include tidal ef¬ 
fects [18-20] also include the merger, and provide better 
accuracy than PN by tuning higher-order vacuum terms 
to numerical relativity (NR) BBH waveforms. Although 
EOB has accurately reproduced waveforms from NR BNS 
simulations [20, 21], here we discuss a new and different 
approach that holds considerable promise for modeling 
tidal interactions during the inspiral. 

The most accurate method of computing waveforms is 
carrying out full NR simulations for BNS and BHNS bi¬ 
naries; see [10, 20-30] for recent work. However, BNS 
and BHNS simulations are computationally challenging, 
since they require solving not only the full Einstein equa¬ 
tions but also relativistic hydrodynamics with a realistic 
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EOS. It is unfeasible to use NR hydrodynamic simula¬ 
tions alone to cover the parameter space given the wide 
range of theoretically possible EOS and NS masses. In 
contrast, BBH systems are easier to simulate with higher 
accuracy. Several large catalogs of BBH simulations and 
resulting waveforms have been compiled [31-37]. 


We introduce here a method we call “PN tidal splic¬ 
ing”, which generates BNS inspiral waveforms from NR 
BBH waveforms by adding tidal interactions derived in 
the PN formalism, effectively replacing the point-particle 
PN terms by the numerical BBH evolution. 


We compare PN tidal splicing to NR using two sim¬ 
ulations generated by SpEC [39], a code developed to 
evolve Einstein’s equations and general relativistic hy¬ 
drodynamics [40, 41]. The first is a new equal-mass 
BNS simulation with 22 orbits before merger [42], and 
the neutron stars were initialized with gravitational 
masses rrii « 1.64M0 and a polytropic EOS with P — 
123.6Mq/ 9^, leading to a tidal deformability of Xi « 
5.7 X 10^®g cm^ s^. The other is an equal-mass, non¬ 
spinning BBH simulation [43] tagged SXS:BBH:0180 in 
the public simulation catalog of the Simulating eXtreme 
Spacetimes Collaboration [35]. Using tidal splicing, we 
add tidal terms to the BBH waveform in an attempt to 
reproduce the BNS waveform. As a test, we also sub¬ 
tract tidal terms from the BNS waveform in an attempt 
to reproduce the BBH waveform. 


Figure 1 shows that the GW phase difference, 
l^'/'Gwjj between the ‘BBH-|-tidal’ waveform and the 
BNS waveform is the same as the difference between the 
‘BNS—tidal’ waveform and the BBH waveform, and both 
are a factor of ^ 3 smaller than the difference between 
the BNS and BBH waveforms throughout the inspiral. 
Thus we can mimic the inspiral of a full BNS simulation 
to within a few tenths of a radian at a fraction of the 
cost. For the BBH waveform, the phase error is estimated 
by the phase difference between the highest two resolu¬ 
tions. The BNS simulation is a combination of spectral 
and finite-volume methods with complicated convergence 
properties; it is unclear how to construct an accurate er¬ 
ror measure [42]. We choose the simple prescription of 
plotting the phase difference between the highest two res¬ 
olutions as a crude error estimate. While the BBH error 
estimate is small, the error estimate in the BNS simula¬ 
tion is as large as the tidal effects themselves. Therefore, 
we cannot yet fully verify the accuracy of tidal splicing 
until more accurate BNS simulations are available. Below 
(cf. Fig. 3) we will estimate the phase accuracy required 
of future BNS simulations for such verihcation. 


II. METHODS 

For nonprecessing binaries, the PN equations for qua¬ 
sicircular orbits read 


dv 

dt 

d(j) 

dt 


=F{v ), 
=v^/M , 


( 1 ) 

( 2 ) 


where F{v) is the ratio of two functions, each known to 
finite PN order in v, and also depends on the binary’s 
intrinsic parameters [44]. Different ways of evaluating 
these equations result in different PN approximants that 
agree to the same PN order in v, but diverge at higher 
orders. We present methods for tidal splicing using two 
different approximants. 

TaylorT^. If F{v) is expanded as a series in v and 
then truncated to the appropriate PN order, then the 
solution is known as the TaylorT4 approximant [45]. For 
TaylorT4, the tidal effects manifest as additional terms 
in the power series for F(v). Equation (1) can be written 


_=F(z;)=Fpp(u) + Ftid(u), (3) 

where Fpp{v) are the point-particle terms, and where the 
additional static tidal terms Ftid(^^) are known to 6PN 
order [15]. 

For inspiraling PN BBHs, F{v) is governed by the 
point-particle terms. PN tidal splicing uses from 

a BBH simulation together with Eqs. (3) and (2) [with 
Ftidiv) set to zero] to compute an accurate version of 
Fpp{v), which we will call Ft^^{v). To do this, we set 
= ^gw/ 2, where (j)GW is the GW phase of the 
£ = m = 2 spherical-harmonic mode of the NR wave¬ 
form. Then Eq. (2) yields 


/ M d(j)GW \ 

V 2 dt J 


( 4 ) 


Given v(t), we compute Fnr(u) = dv/dt using finite dif¬ 
ferencing. Assuming v{t) is monotonic, we can write 
Fnr(u) as a single-valued function of v. 

Using this FNR(r’) in place of Fpp{v) in Eq. (3), we 
then re-solve Eqs. (3) and (2), including the tidal terms 
Atid(^')) to generate a waveform for a binary that includes 
tidal interactions. We express the orbital evolution of the 
new binary in terms of a new time coordinate i. From 
the analytic expression for Ftid{v) [15] and Eq. (3) we 
write a differential equation for t: 


dv Enr(u)-I-F tid(u) ■ 

Integrating this expression and inverting yields the func¬ 
tion v{t) corresponding to the spliced waveform. 
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FIG. 1. Phase difference between gravitational waveforms as a function of time, for an equal-mass binary of nonspinning 
compact objects. Differences are shown between BNS and BBH waveforms (black), between a BBH waveform with TaylorT4 
tidal terms added and a BNS waveform (blue), and between a BNS waveform with TaylorT4 tidal terms subtracted and a BBH 
waveform (red). The red and blue curves nearly coincide. Also shown are the phase differences between BBH and point-particle 
TaylorT4 waveforms (solid magenta) and between BNS and tidal TaylorT4 waveforms (dashed magenta). The numerical error 
in the BBH waveform (dashed red) and an estimate of the error in the BNS waveform (dashed blue) are also shown. All 
waveforms are aligned with the BNS waveform according to [38]; the alignment time window encompasses a 5% change around 
a GW frequency of 280 Hz for a total mass of M = 2 x l.OdM©. The blue and red curves are smaller than the black curve by a 
factor of ~ 3, demonstrating that tidal splicing can generate a BNS waveform from a BBH waveform and vice versa. The large 
error in the BNS waveform prevents us from fully measuring the accuracy of tidal splicing. 


The phase of the spliced waveform, (/>Gw(^) is com¬ 
puted by integrating Eq. (2): 

bGw(i) = ^^ vit}^dt, (6) 

where tmm is the start of the numerical waveform. 

In TaylorT4, the amplitude of the waveform is a func¬ 
tion of V only, with no explicit time dependence [46]. 
So here we assume that the amplitude of the original 
NR waveform Anr)!) is likewise a function of v only, 
so that we can write Anr)^) = ANR(t(f)). We then 
use u(t) to express this amplitude in terms of i. In 
other words, the amplitude of the resulting waveform is 
m = . We generate a BBH waveform from 

a BNS waveform by the same method, except we subtract 
instead of add Etid('c) in the denominator of Eq. (5). 

We require v{t) to be monotonic so that F{v) is sin¬ 
gle valued. To remove high-frequency numerical noise, 
the derivative in Eq. (4) is computed with a third order 
Savitzky-Golay filter [47] with a window size of « 48.5 /xs. 
This is sufficient when adding tidal terms to the BBH 
waveform considered here. However, when testing our 
method by subtracting tidal terms from a BNS waveform, 
the phase of the BNS waveform considered here [42] has 
large enough oscillations in v{t) that stronger smooth¬ 
ing is needed. We proceed by first subtracting the phase 


of the TaylorT4 waveform from that of the BNS wave¬ 
form, expanding this difference in Chebyshev polynomi¬ 
als, truncating the Chebyshev expansion to n = 35, and 
adding back the phase of the TaylorT4 waveform. We 
find that the difference between the smoothed and un¬ 
smoothed phase of the BNS waveform is less than 3x 10“^ 
radians. 

As discussed above. Figure 1 displays the phase dif¬ 
ferences between NR and tidally spliced TaylorT4 wave¬ 
forms. We now examine how well pure PN waveforms 
agree with NR waveforms. The magenta solid and dashed 
curves in Fig. 1 show phase differences between TaylorT4 
and BNS or BBH waveforms. The point-particle Tay- 
lorT4 waveform does an excellent job of reproducing the 
phase evolution of the BBH waveform, about at the level 
of the BBH numerical error. However, while TaylorT4 
is surprisingly accurate in the inspiral for equal-mass, 
nonspinning systems [45, 48], this does not hold true in 
general [49-51]. Tidal splicing should be applicable to an 
arbitrary BNS/BHNS system with spins and/or unequal 
masses, where there may not be an accurate PN approx- 
imant. References [11, 12] showed that uncertainties in 
the PN waveforms are one of the largest sources of error 
for tidal parameter estimation, and conclude that more 
accurate waveforms are needed. 

TaylorF2. If Eqs. (1) and (2) are instead converted to 
the frequency domain (ED) using the stationary phase 





4 


approximation before expanding the series, the approx- 
imant is called TaylorF2 [52]. TaylorF2 waveforms are 
expressed in the FD, and can be written 

^(/) = ^(/)exp , (7) 


where A{f) is real and 4'(/) is the Fourier phase as a 
function of the GW frequency / = v^/{ttM). For point 
particles, ^(/) = ^pp(/) is known for nonspinning sys¬ 
tems to 3.5PN order [52, 53]. For tidally deformable ob¬ 
jects, we write 4'(/) = 4'pp(/) -f ^tid(/), where ^tid(/) 
has been calculated up to 7.5PN order, with the excep¬ 
tion of a few unknown constants [7, 19]. Here we include 
both 6PN tidal effects and 7.5PN tidal effects, setting the 
unknown constants to 0 as was done in [13]. 

To add the static tidal terms to an existing BBH wave¬ 
form, first the Fourier transform of the waveform /inr(/) 
is computed. The early portion of the waveform is win¬ 
dowed using a Planck taper [54] while the merger and 
ringdown provide a natural windowing for the late por¬ 
tion. We then compute ^nr(/) and Hnr(/) by decom¬ 
posing according to Eq. (7). The spliced Fourier phase 
is then 4'(/) = ^nr(/) + ^tid(/)- Because the known 
tidal terms do not affect the amplitude Hnr(/), the new 
waveform is then 


Hf) = -4nr(/) exp 4 'nr(/) + ^tid(/) ) • (8) 


No smoothing of the numerical waveforms is needed for 
TaylorF2 splicing, unlike the TaylorT4 case. 

Since the PN approximation breaks down for high fre¬ 
quencies, we impose a high frequency cutoff which we 
choose to be /isco = l/(6^/^7rM) = 1338Hz, the GW 
frequency corresponding to the innermost stable circu¬ 
lar orbit of a Schwarzschild black hole of mass equal to 
the total mass of the system. It has been shown that 
for BNS systems, /isco is approximately the merger fre¬ 
quency [55]. The starting frequency of the NR BNS wave¬ 
form after windowing is ^ 285 Hz. 

We estimate the error of the spliced waveforms by an¬ 
alyzing the phase differences in the time domain after 
taking the inverse Fourier transform. To avoid jump dis¬ 
continuities in the Fourier phase, we roll off the effect of 
^tid(/) from /isco to 2 X /isco with a cosine window. 
While this will contaminate the higher frequency content, 
this should allow the lower frequencies of the inspiral to 
be mostly unaffected. After the inverse Fourier trans¬ 
form, the time domain waveforms are aligned in a 10% 
region around 300 Hz. The phase differences are shown 
in Fig. 2 and are similar to Fig. 1. With the exception 
of the last ^ 3 ms of the waveforms, which are affected 
by the high frequency contamination, all of the spliced 
waveforms maintain phase differences under 0.1 radians 
during most of the inspiral, below the difference between 
the BNS and BBH waveforms. It is not clear why the 
6PN terms approximate the tidal effects better than the 



Time (s) 

FIG. 2. The phase difference \S(j>aw{t)\ as a function of 
time for waveforms spliced with TaylorF2. Differences are 
shown between a BNS and a BBH waveform (black), between 
a BBH+tidal and a BNS waveform (blue), and between a 
BNS—tidal waveform and a BBH waveform (red) at the 6PN 
(solid) and 7.5PN (dot-dashed) orders. Only the time after 
the windowing function is shown here, resulting in a shorter 
time axis here than in Fig. 1. The late-time noise is an ar¬ 
tifact caused by inverse Fourier transforming the unphysical 
high-frequency behavior of ^tid{f)- At both PN orders, tidal 
splicing can generate a BNS waveform from a BBH waveform 
and vice versa. 


7.5PN terms; it may be because we zero the unknown 
constants in the 7.5PN expression. 


III. DISCUSSION 

We have shown that PN tidal splicing of BBH wave¬ 
forms can produce inspiral waveforms for nonspinning 
BNS systems. This method should easily generalize to 
objects with spins and to BHNS systems. Once a BBH 
waveform is generated for a particular mass ratio and spin 
configuration, it should be easy to produce BNS/BHNS 
waveforms via PN tidal splicing for any EOS simply 
by adjusting the tidal parameters Xi, allowing the en¬ 
tire tidal parameter space for inspiral waveforms to be 
spanned. 

The accuracy of this method is limited by that of the 
PN tidal terms. In particular, additional finite size ef¬ 
fects not captured by the current static tidal PN terms 
can influence waveform amplitude and phase, and dy¬ 
namical tidal effects can also contribute to the phase 
evolution [56], especially as the NSs approach merger. 
This method in principle can be improved with better 
PN tidal terms. Unfortunately, it is currently difficult to 
fully measure the accuracy of tidal splicing until higher- 
accuracy many-orbit BNS simulations are available for 
multiple masses and EOS. 

Figure 3 estimates the accuracy needed for equal mass. 
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FIG. 3. Phase difference between equal-mass, nonspin¬ 
ning BBH and ‘BBH-|-tidar waveforms. Each horizontal slice 
through this plot shows the phase difference as a function of 
time for a particular dimensionless deformability Ai/mf. For 
our BNS simulation, Ai/mf « 453. Contours show selected 
values of the phase difference. A BNS simulation starting at 
dimensionless time t/M ~ —4500 would need phase errors 
smaller than the values shown here in order to measure tidal 
effects. Even more accurate BNS simulations would be needed 
to measure the accuracy of the tidal splicing method. 
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nonspinning BNS simulations to see the tidal effects on 
the inspiral phase of the waveform. Even smaller BNS er¬ 
rors would be necessary to constrain the accuracy of tidal 
splicing. We chose the start time in Fig. 3 so that the 
inspiral spans a large enough frequency range for aLIGO 
to recover 97% of the information about A^, according 
to the analysis presented in Fig. 3 of [7]. We assume 
M = 2.8Mq (corresponding to a prototypical NS mass 
of I. 4 M 0 ) and an upper frequency cutoff of /isco- 

An alternative to computing tidal terms to a higher PN 
order is to resum them in some way, as is done in [20, 21, 
57] in the context of EOB. It is not clear how to do this 
with tidal splicing. 

Additionally, the merger/ringdown cannot be mod¬ 
eled with splicing alone, especially for BNS mergers and 
BHNS systems that undergo tidal disruption. One pos¬ 
sibility is to combine an analytic waveform in the very 
early inspiral with a spliced BBH waveform in the mid to 
late inspiral and then with a waveform from a full hydro- 
dynamical simulation for the merger and ringdown, to 
create a “tribridized” waveform. This might reduce the 
need for expensive hydrodynamical simulations lasting 
many orbits. If necessary, surrogate models [43, 58, 59] 
that cover the parameter space including the EOS can 
be forged from spliced BBH waveforms. 
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